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0^ ■ ABSTRACT 

Distribution functions (DFs) for dynamically warm thin stellar disks residing in arbi- 
P5 ' trary axisymmetric potentials are presented which approximately reproduce pre-described 

^ ■ surface-density and velocity-dispersion profiles. The functional form of the DFs is obtained 

by 'warming-up' a model made entirely of circular orbits. This can be done in various 
ways giving different functional forms for the DP. In the best case, the DF reproduces 
the pre-described profiles to within a few per cent for a typical case reminiscent to the 
old stellar disk in the Milky Way. This match may be improved to about one per cent or 
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00 '. better by a simple iterative method. An algorithm is given to draw phase-space points 

^ \ randomly from the DFs for the purpose of, e.g., A^-body simulations. All the relevant 

Q ■ computer programs are available from the author. 
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1. Introduction 

All the dynamical properties of a collisionless stellar sys- 
tem are prescribed by its phase-space distribution func- 
tion (DF) in conjunction with the underlying gravita- 
tional potential $. The knowledge of a simple DF is 
often vital for studies of the dynamical state of stellar 
systems, their equilibrium properties and time-evolution. 
In this paper, simple DFs for thin axisymmetric disks are 
given. These DFs may be used in studies of, e.g., spiral 
structure, the effects of non-axisymmetric perturbation 
on the velocity distribution in the solar neighborhood, 
and to seed initial conditions for A^-body simulations. 

For a thin axisymmetric disk, ^ — ^{R) and the col- 
lisionless Boltzmann equation tells us that in dynamical 
equilibrium the stellar DF must be a function of the stel- 
lar energy and angular momentum alone, f — f{E,L). 
Thus, compared to the more general case of three- 
dimensional axisymmetric or non-axisymmetric systems, 
where the DF may also depend on non-classical integrals 
of motion, the DFs of thin axisymmetric disks are simple 
and allow for mathematical manipulations. It is possible 
to give a DF /(-B, L) which implements some pre-defined 
properties, for example, a certain surface density pro- 
file S(i?), in a given underlying gravitational potential 
<I>(i?)Q Unfortunately however, this DF is highly non- 
unique: every one of the infinitely many ways to expand 
E(i?) as a function of both R and <&(i?), and thus to 
obtain a so-called surface-density partition I](i?, $), de- 



termines uniquely the part of the DF even in L (Kalnajs 
1976). This even part of the DF, /+, ca n be computed 



(usua lly numerically) by a path-integral ( Hunter fc Qian 
1993 ). Additional constraints, for example, a certain 



radial- velocity-dispersion profile crfl(_R), restrict the pos- 
sible surface-density partitions and thus DFs, but do not 
remove the general degeneracy (obviously, the DF being 
a function of two variables cannot be uniquely deter- 
mined by a finite set of functions of one variable) . 

Given this degeneracy among possible disk DFs, there 
are two different ways to construct specific dynamical 
models. The first approach follows the above lines in 
specifying a particular surface-density partition and then 
solving for the implied f+{E,L). This technique has 
the advantage that the pre-defined properties are exactly 
satisfied by the resulting DF, but the disadvantage that 
the DF can rarely be given in closed form, which tends 
to be complicated. 

The second way to a disk DF is to specify a particular 
functional form, that results only in approximate agree- 
ment with the pre-defined properties. This approach has 



^ In this paper, we will not restrict ourselves to self-consistent mod- 
els in which the gravitational potential is itself generated by the 
surface density given through the DF, but allow for a general #, 
which may partly be due to the stellar disk. 



the advantage that the DF is always of a simple func- 
tional form, which often is of great importance in dy- 
namical modeling. Another bonus is that the degeneracy 
among DFs with given surface density and velocity dis- 
persion is removed by choosing a particular functional 
form for the DF itself rather than the partition func- 
tion. This enables one to choose the functional form for 
f{E, L) on the basis of astro-physical arguments. On the 
other hand, in the first approach the degeneracy is often 
exploited by choosing the surface-density partition such 
that the resulting DF is analytic. 

The aim of this paper is to improve and extend on 
what has been done before on that second approach to 
warm-disk DFs. In Section 0, warm-disk DFs are con- 
structed by warming-up the DF for a completely cold 
disk in which all stars on circular orbits. In Section y, the 
ability of four types of warm-disk DFs to reproduce the 
predefined surface-density and velocity-dispersion profile 
is assessed and a simple algorithm is given to improve on 
this ability. Section [J compares the velocity moments up 
to fourth-order for the four DFs introduced in Section ||, 
while Section || gives an algorithm for drawing from the 
DFs samples of phase-space points such as might be used 
for numerical simulations. 

2. Recipes for Warm-Disk Distribution 
Functions 

We seek DFs f{E, L) which approximately create given 
targets for the surface density S(-R) and radial velocity 
dispersion aR{R) in a given potential $(i?). The basic 
technique is to modify the DF for a completely cold disk. 

2.1. The Distribution function for the cold disk 

Let Ec (L) be the energy of the circular orbit with angular 
momentum L (see, e.g.. Appendix A of Dehnen 1999, 
hereafter paper I, for details of the properties of circular 
orbits), then 



fcoM{E,L) = U,{L) S{E - E,{L)), 



(1) 



describes a DF which allows for circular orbits only. 
From 



M^{2^f j d'j f, 



(2) 



which may be derived using d^a; d\; = d^J d^ with J = 
{Jr, L) denoting the action variables, and {dE/dJn)L = 
cjfl,, we have 



M = (27r)^ dE dL f{E, L) uor^{E, L), 
and thus 

dM = (2^)Veirc(i)/«(i) dL, 



(3) 



(4) 
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where I have used the fact that on circular orbits the 
radial frequency lj^ equals the epicycle frequency k. On 
the other hand we have 



dM = 27ri?S(i?)di?, 



(5) 



and equating (Q) with (||) gives, with dLc/dR — Rn/^ 
(cf. equations A2 to A6 of paper I), fdidL) = /i (i?c(-^)) 
with 



h{R) 



2n 



(6) 



Here, Rc{L) is the radius of the circular orbit with an- 
gular momentum L, while 7 = 2n/K with the circular 
frequency fl{R). 

2.2. The Shu distribution function 

In order to allow for non-circular motions, the (5- function 
in the cold-disk DF (|l|) must be replaced by some func- 
tion of finite width. Since the epicycle energy E — Ec{L) 
is non-negative, a simple choice is an exponential corre- 
sponding to a Maxwellian velocity distribution 



5{E-E,{L)) 



1 



■ exp 



Ec{L) - E 



(7) 



Using i?L =Rc{L) also as argument of cr|j, this leads to 
the well-known DF (Shu 1969) 



fshn{E,L) 



j{RL)^iRL) 

2TTa%{RL) 



exp 



Ec{L)-E 

<^URl) 



(8) 



2.3. New warm-disk distribution functions 

Instead of (Q) , we could have started with the alternative 
ansatz for a cold-disk DF: 



fcoid{E,L) = /eirc(^) s(n{E)[L,{E) - L] 



(9) 



With dE,/dR = RKy2, one finds UrdE) = h{Rc{E)), 
where Re (E) is the radius of the circular orbit at energy 
E. As /shu was derived from O), we may derive a new 
warm-disk DF from (pi): 



Jncw\E, L) 



-f{RE)J:{RE) 

2ti a\{RE) 



exp 



n{RE)[L-LdE)] 
<{Re) 



(10) 



with Re = Rc{E). There is actually large freedom in de- 
signing more warm-disk DFs by combining the various 
possibilities for the arguments of /i, cr|j, and the expo- 
nential. Two other possible DFs are 



UE.L) 
ME,L) 



-f{RE)nRE) 
2TTa%{RE) 

jiRL)nRL) 
27ra|(i?i) 



exp 



exp 



EciL)-E 



aliRE) 
n{RL)[L-L,{E)] 



^URl) 



(11) 

(12) 



2.4. Physical Motivation 

Shu (1969) motivated his DF by the idea that an ini- 
tial process of violent relaxation formed the distribution 
function of stellar disks. He also assumed that the ini- 
tial distribution of angular momenta is not changed by 
this process. However, we now know, that these ideas 
are incorrect. Disk stars gain their random motions by 
scattering off molecular clouds, spiral waves, satellite 
galaxies and other agents. Therefore, the velocity dis- 
persion grows continuously, which is also reflected in the 
observed increase of random motions with stellar age in 
the solar neighborhood. Such scattering processes do not 
change the value of a star's Hamiltonian in the frame co- 
rotating with the scattering agent. This Hamiltonian is 
equal to the Jacobi energy Ej = E — loL with to the an- 
gular frequency of the scattering agent. Consequently, 
both the stellar energies and angular momenta are af- 
fected by disk heating. Thus, we may describe the DF 
as a warmed-up version of a DF made from new-born 
stars. Stars are born already with some finite velocity 
dispersion, which is small compared to that gained sub- 
sequently, and we may safely assume the DF of new-born 
stars is dynamically completely cold. The warming-up 
may be described by an exponential in the radial action 



Jji (Binney 1987) 



/ « fi{Ri) 



i(Ri) 



exp 



^rJr 



<(.Ri) 



(13) 



Here, i?i = i?i {E, L) is some measure of the mean orbital 
radius that ensures the prescribed surface density to be 
closely matched by the DF. The radial action may be 
estimated from the classical epicycle theory to be 

ujRjR^E-EdL), 

but, as demonstrated in paper I, 

ujRjR^n{E)[LdE)-\L\] 

is much better an approximation for near-flat rotation 
curves. Using either of these approximations for ujuJr 
and either Rl or Re for i?i, any of the four DFs pre- 
sented above can be derived from equation (|l3|). 

There are two arguments in favor of /now being the 
most useful of these four DFs. First, the mean orbital 
radius i?(-E, L) is much better approximated by Re than 
by Re (in general: Rl <Re^R)- Similarly, k{Re) ap- 
proximates the radial frequency lor, which we replaced 
by K in the derivation of /i, far better than k{Rl), see, 
e.g. Fig. 8 of paper I. Therefore, one expects DFs with 
Rl ~ Re to create surface density and velocity disper- 
sion in better agreement with the target functions than 
DFs with Ri=Rl. Second, 0(L — Lc) as argument of 
the exponential naturally extends to negative L, while an 
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Fig. 1. — Logarithmic deviations between the target surface density and velocity dispersion (S, aji) and those actually created by the 
warm-disk DF (S/, <^R,f) in power-law potentials with circular speed Vc cc R^ . Solid lines are for /shu, dashed for /now, dot-dash for /a, 
and long-dashed for /b- In all four sub-figures, both S and aji decay exponentially with scale lengths R^ =3Rb- In (a), (b), and (c), the 
central velocity dispersion is O.bvg, while it is twice as large for (d). The rotation curve is flat in (a) and (d), slightly falling in (b) and 
slightly rising in (c). Notice the different scales of the y-axes. 



exponential in E — Ec does not allow for this possibility. 
This is significant insofar, as the tail of stars with small 
negative L is known to be important for supporting disks 
against dynamical instabilities. 



3. Assessing the Distribution Functions 

Allowing for non-circular motions, i.e. replacing the S- 
function in (nl) or M) with some function of finite width, 
leads at any position i? to a mixture of stars originating 
from different radii Ri. As a consequence, the surface 
density S/(i?) that is actually created by the DF devi- 
ates from the target I]{R). Similarly, the radial velocity 
dispersion aRj{R) due to the DF is expected to differ 
from the target afi{R). On the one hand a simple and 
generally applicable functional form for the DF is very 
desirable, on the other hand it is equally important for 
the DF to re-produce these targets as closely as possible. 
In j j3.l| , we assess the ability of the DFs introduced in 
§0 to satisfy this latter goal, while j |3.2| presents a sim- 
ple numerical method by which the match to the target 
functions can be vastly improved. 



3.1. How Well do the DFs Re-Produce their 
Targets? 

We now specify the target functions T,{R) and (Jr{R) and 
compare these with the numerically evaluated velocity 
moments resulting from the various forms of warm-disk 
DFs. We restrict ourselves to exponential surface-density 
and velocity-dispersion profiles 



aR{R) = ao e~^/^' 



(14) 
(15) 

with scale radii Rg and Ra- The gravitational poten- 
tial is assumed to be a power-law with circular speed 
curve Vc — vo{R/RsY (cf- Appendix B of paper I). Fur- 
thermore, we only consider models with R„ = 3i?s re- 
sulting in cr^ ex E corresponding roughly to the situation 
in the Milky- Way disk. For flat [(3 — Q), slightly rising 
(/3 = 0.2) or falling {13= — 0.2) rotation curves, and for 
the four warm-disk DFs introduced in the last section. 
Figure plots the logarithmic deviations of surface den- 
sity and velocity dispersion for ctq = O.Suq (Fig. 0a-c), cor- 
responding to cr(i?o ~2.5i?s) «0.2uo reminiscent of the 
old-stellar disk of the Milky Way, and doubly as much, 
ao=woXor/3 = Oonfy(Fig. Id). 

At R^Rs in all four cases, Shu's DF shows the 
strongest deviation both from the target surface density 
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Fig. 2. — The same as in Fig. g^, but for the DFs fshu (a) 
and /new (b) after none (thin), one (dotted), two (dashed), three 
f lo ng-das hed') . and four (bold dot-dash) iterations of the scheme 
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and velocity dispersion, while /now gives the best match. 
On the other hand, at small radii /new, and to a smaller 
extent /a, creates strongly rising E/ and ctrj. The likely 
cause of this effect are eccentric orbits contributing sig- 
nificantly to the surface density near their peri-centers, 
i.e. at small radii. In most applications, however, this 
should be no problem, since (i) there is negligibly little 
mass in these artificial cusps, and (ii) the central parts of 
disk galaxies are dominated by dyna mical ly hot bulges. 
Moreover, the technique described in §B^ below may be 
used to suppress these artificial cusps. 

The case with an target of ctq = vq (Fig. |I|d) clearly 
represents more a hot than a cold disk {aj^ + a'^>v'^ out 
to i?~ I.3i?s for /new), while the DFs were only meant 
to work for moderately warm disks (i.e. ordered motions 
dominating). Nonetheless, the moments created by /new 
match the target functions reasonably well. The other 
two new warm-disk DFs, /a and /b, are less useful than 

/new- 

From Figure n^, we see that the DFs are much better 
in re-producing their target functions for rising than for 
falling rotation curves. Presumably, this is because for 
a fixed (jr[R)^ the ratio ur/vc is, at most radii, smaller 
for rising than for falling rotation curves. A smaller ratio 
^r/Vc implies that the disk is colder and the DFs, which 




Fig. 3. — As Fig. Q but for the stellar disk and rotation curve 
as in Fig. lid (which is hot inside ~ l.S-Rj;) and after none (thin), 
one (dotted), two (dashed), five (long-dashed), and seven (bold 



2.1, Notice the different 



dot-dash) iterations of the scheme of 
scales on the j/-axes. 



are exact only in the cold limits, give a better match. 

3.2. Improving the distribution functions 

The results of the last subsection show that DFs with 
(jR / only approximately re-produce their target sur- 
face density and velocity dispersion. The goal of this 
section is to get a much better match with some given 
surface density 5](i?) and velocity dispersion cr_R(i?), by 
distinguishing between the targets and the parameter 
functions that enter the functional form of the DF. That 
is, we replace the functions S and gr, in the defini- 
tions of the warm-disk DFs with S' and a^, which are 
to be constructed such that the resulting moments of 
the DF closely match the target: Ey « E ^ E' and 
ctrj^ctr ^ o-^. 

3.2.1. A Simple Algorithm to Construct E' and a'j^ 

The simplest technique to optimize E' and cr^ is to start 
with E' = E and cr^ ~ a and iterate the following steps. 

1. Compute Ey and cfrj due to the DF. 

2. Multiply E' with E/Ey and cr^ with CR/aRj. 

Figure shows the logarithmic deviations between target 
and actual E and aR that have been obtained after up 
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to four iterations for the same stellar disk and potential 
as in Fig. nk. There is clearly a dramatic improvement 
already after one iteration: the average relative deviation 
has dropped by up to a factor ten. Figure || is similar to 
Fig. |[ except that the target a^ is twice as large (the 
model used in Fig. ||d) such that random motions are 
actually dominant inside about 1.3i?i;. This time, the 
improvements are modest, and only in the case of /now 
do the iterations reduce the deviation from the targets 
to a few per cent. 

It seems, that after a few iterations no further im- 
provement is achieved, i.e. the algorithm does converge 
to a model that definitely deviates from the target func- 
tions. For modest velocity dispersions this deviation is 
less than about one per cent, but up to a few per cent or 
even more, depending on the DF, for stellar disks with 
significant random motions. From our explanation for 
the algorithm below, this is not astonishing, since the 
algorithm is very crude indeed; on the contrary, it is 
suprising that it works so well. 

3.2.2. Why Does This Algorithm Work? 

The simple algorithm given above can be derived as fol- 
lows. The integrations over velocity space that yield the 
moments may be cast into the form 



E/(i?) = /"s'(i?i) Xs(i?|i?i)di?i, 
f^fiR)= f ii'{Ri)K^{R\Ri)dRi 



(16a) 
(16b) 



with n = u^E. The kernels Ky: and K^^ are implicitly de- 
fined by the form of the DF, but cannot be obtained ex- 
plicitly. However, epicycle theory tells us that these Ker- 
nels should have a typical width of order a'j^{Ri) / k{Ri), 
which in turn is of the order {a'j^/vc)RiRi. Equation ( |l6[) 
defines an inverse problem for finding the primed quanti- 
ties given the unprimed ones. Clearly, to yield a physical 
meaningful DF, the primed quantites must not become 
negative, and a simple technique to solve for E' and tr^ 
whilst ensuring positive definiteness is the Richardson- 
Lucy algorithm. This consists of iterating on (|l^) and 
the correction step 

E'(i?i) ^ E'(i?i) X f ^^ Ky{R\Ri) dR, (17a) 

fi'{Ri) -> fi'{Ri) X J 4^ K^{R\Ri)dR. (17b) 

The algorithm given above emerges when replacing the 
kernels in (^TI) with J-functions. Clearly, this is a very 
crude simplification, but computing the integrals in (QTh 
is non-trivial, and furthermore, the scheme works rea- 
sonably well, as demonstrated above. 



4. The Kinematics 

As mentioned in the introduction, the distribution func- 
tion f{E,L) for a stellar disk is not uniquely specified 
by its surface density, velocity dispersion, and the grav- 
itational potential. Rather, there is an infinite number 
of possible distribution functions with identical E^ and 
o'Rj, but with differences in, e.g., other velocity mo- 
ments. Here, we will study the kinematics, i.e. the ve- 
locity moments up to order four, that emerge from the 
various warm-disk DFs when applied to the same stellar 
disk. 

For the case of an exponential disk with i?(j = 3i?s 
and cro = 0.5uo (as in Fig. ||a-c). Figure ^ plots various 
kinematical quantities resulting from the four different 
DFs of §||. Figure ^ differs only in that the DFs are 
those obtained after four iterations of the algorithm of 
§3.2.1. Obviously, in some cases this rather crude algo- 



rithm tends to create DFs with somewhat "ragged" kine- 
matics. 

The top panels show the run of the asymmetric drift 
velocity, which usually peaks between 1 and 2 iJs for 
these models. The iterated DFs show almost identical 
behavior in their asymmetric drifts. 

The second panels from top show the ratio between 
the kinetic energies in azimuthal and radial random mo- 
tions, which in the solar neighborhood is measured to be 
about 0.42 (cf. Dehnen & Binney 1998). In the limit 
of small velocity dispersions, we know from Oort's work 
that 



lim 



-0 at 



1 



d In Vr 



d\nR 



-B 



A^B 



7 



(18) 



where A and B are Oort's constants. Indeed, at large R, 
this value is reached by the moments of the DFs. How- 
ever, at radii R^ 2i?s, all four DFs predict cr'i/a'j^ to be 
significantly larger than this limit. This is a known result 
for exponential disksfl (cf. Evans fc CoUett 1993| , [Cuddc- 
ford fc Binney 1994|). As the third panels from top of 



Figures ^ and H show, the v^ distributions are always sig- 
nicantly skewed in the sense that they have a low-w^ tail. 
This tail is caused by stars originating from Ri < R with 
large eccentricities - because of the exponential E(_R) and 
aii{R) relations, there are many fewer stars with large 
eccentricity originating from i?i > R, which would create 
a high-Vcj, tail. Such a skewness in the azimuthal motions 
is indeed observed in the solar neighborhood, cf. Dehnen 
(1998). 

It is only in the kurtoses, i.e. moments of the fourth 
order, that the four DFs appear to differ significantly. 



^ This fact is not very well known, though, and it cannot be stressed 
enough that Oort's equation (Hh) is correct only in that limit and 
must not be applied to the situation in the solar neighborhood, 
even though this has been done many times in the past. 
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Fig. 4. — Kinematics vs. R for the warm-disk DFs /shu (solid), /new (dashed), /a (dot-dashed), and /b (long dashed) for an exponential 
disk with iJo-= SiJj; and (To = 0.5i)o in a slightly falling (/3 = — 0.2), flat (/3 = 0), and slightly rising (fi = 0.2) rotation curve (same as in 
Fig. Ba,b,c). 



While the kurtosis for the azimuthal motion is mostly 
positive, that of the radial motion may have either sign. 

5. Sampling of Phase-Space Points 

For A^-body simulations and similar studies, it is of im- 
portance to draw phase-space points (i?, i?, 0, 0)^, i = 
0, . . . iV from the distribution function. One is often in- 
terested not in randomly distributed phase-space points 
but in a more regular distribution, where the points are 
equidistant in the phases (angle variables). Such a distri- 
bution has minimal noise, a fact important in numerical 
studies of stability properties ( "quiet start" technique) . 
From the general form (O) of the warm-disk DFs pre- 
sented, it is straightforward to derive a procedure for 
sampling [E^ L) as follows. First, determine i?i accord- 
ing to E(i?), and translate it to either E = Ec{Ri) or 
L = Lc{Ri) according to the functional form of the DF. 
Second, sample the exponential in radial action and com- 
pute the remaining integral: E or L. Here is a general 
algorithm that allows for parameter functions E' and a' 
different from the targets and accurately corresponds to 
the DFs of §|. 



1. Choose a radius R randomly from the cumulative dis- 
tribution 



P{R) = f R' S(i?') dR' 
Jo 



(19) 



and determine L = Lc{R) for fshu and fb or E = Ec{R) 

for /now and /a. Evaluate the correction factor gi — 

S'(i?)/S](i?). 

In the case of /b, the sign of L must be chosen to be 

positive or negative with relative probabilities 

l:exp[~2f](i?i)|i|/a|j(i?i)]. 



(20) 



2. Choose ^ G (0,1) randomly and determine 

E = E,{L)~a'l{R)\n^ for /ghu, 

L = L,{E) + a'^{R) lnC/n{R) for /„ew, 

L^L,{E + a'^iR)ln^) for /a, 

E = E,{\L\-a'^{R)\ii^/n{R)) for /t. 

HE ^ [$(0),$(oo)] or L ^ \~L,(E),L,(E% go back 
to step 1 for another try. 

3. Integrate the orbit with these values of (i?, V) over one 
radial period Tr, compute ujr(E, L) = 27r/Tij, and 
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Fig. 5. — Same as Figure H but for DFs after four iterations of the algogithm of §3.2.1 

evaluate the correction factor 172 — K{R)/oJii{E, L). 



4. Sample Nsam phase-space points randomly from <j) G 
[0, 2tt) and t G [0, Tr), use a table made during orbit 
integration to obtain R{t) and R{t); (p^L/R"^. Here 
iVsam is either of the two integers next to gi g2 No^b, 
chosen with probabilities such that the mean is equal 
to gig2Norh- 

5. Iterate steps 1 to 4 until the desired number N of 
points is sampled. 

Here, iVorb is, approximately, the number of points cho- 
sen per orbit, and is usually between 1 and N^/'^. If 
noise is to be minimized, e.g. if initial conditions for a 
quiet start are needed, use A^orb 3> 1 and quasi-random 
numbers or even equidistant points, otherwise pseudo- 
random numbers. If only (i?, L) shall be sampled, set 
-^orb = 1- Note that for /sim and /b the correction factor 
g2 can become quite large for eccentric orbits, since these 
DFs are based on estimating cj^ by k(£), which goes to 
infinity as i — > in centrally non-harmonic potentials, 
as, e.g., in the Milky Way. In the case of /b, one first 
determines \L\ and has to choose the sign afterwards. 

I have checked the correctness of the above algorithm 



numerically by computing Ey and anf from 10^ sampled 
phase-space points: they agreed within their statistical 
uncertainties with the results from direct numerical in- 
tegration. 

6. Conclusion 



6.1. Summary 

Distribution functions for dynamically warm stellar disks 
may be considered as warmed-up versions of completely 
cold disks, in which all stars rotate on circular orbits. 
In fact, this picture for the DF of a stellar disk may 
be considered a good description for its formation his- 
tory. In Section 0, I presented four different functional 
forms for warm-disk DFs f{E^ L), given the gravitational 
potential $(i?) and prescriptions for the surface-density 
and velocity-dispersion profiles, S(i?) and aii{R). One 
of them, /shu, was introduced by Shu in 1969 and has 
already been used in studies, e.g. of the solar neighbor- 
hood kinematics (see Bienayme 1999 for a recent exam- 
ple). The other three DFs have not been described in 
the literature so far. 



Distribution functions for disks 



In Section ^ the surface-density and velocity- disper- 
sion profiles produced by the DFs as moments were com- 
pared to the target profiles. The comparison shows that, 
for a typical case reminiscent to the old stellar disk in the 
Milky Way, Shu's DF, /shu results in the largest devia- 
tions between moments and targets, while the new DF 
(10) gives deviations that are about three times smaller. 
This new DF even gives an useful model when the stel- 
lar disk is dynamically hot within about one scale radius. 
Actually, the deviation of the moments from the target 
profiles S(i?) and (jji{R) may be considerably reduced 



by lalluwiiig the currespuiidliig parauicLer fuiicLiuiis ap- 
pearnig iii the defliiiLiuii of the DFs Lo dilTei' from these 



targets. In Section 3_^, a simple algorithm for this pur- 
pose is presented, which after only four iteration reduces 
the deviations in S and an. to well below 1 per cent for 
a typical case. Again, the DF /now is best suited for this 
game, as it leads to the smallest deviations also in the 
iterated case. Another advantage of /new is that is nat- 
urally extends to negative L, while /shu make the rather 
unnatural assumption that there are no stars with i < 0. 

In Section 0, I have compared the velocity moments 
up to fourth order of the four DFs iterated and non- 
iterated. For all four DFs, the ratio of azimuthal to radial 
velocity dispersion squared, cr^/cr^ is significantly larger 
than i[l -f dlntic/dlni?] for R larger than about 2 disk 
scale radii. The value_i[l -|-dlnUc/dln_R] is expected 
from Oort's equation dlq), which in turn is based on the 
epicycle approximation whereby ignoring the strong gra- 
dients in an exponential disk. In the solar neighborhood, 
(T^/cr^ ~ 0.42 is observed, i.e. smaller than what we ex- 
pect for a near-flat rotation curve. This deviation must 
be caused by a wrong assumption in the models, most 
likely that of axisymmetry - an investigation in terms of 
the influence of the Galactic bar is the subject of a future 
paper. 

Finally, in Section ^ I give an algorithm for the sam- 
pling of phase-space points {R,pfi,4), L) from any one 
of the four DFs. This should be particularly useful for 
the creation of initial conditions for A^-body simulations, 
either of the evolution of disk galaxies (bar formation, 
warps, etc.) or of mergers involving disk galaxies. 

6.2. Thick Disks 

An obvious extension of the thin-disk DFs presented in 
this paper are DFs describing stellar disks with some fl- 
nite thickness. For this purpose, one needs a handle on 
the vertical motion of the stellar orbits. Unfortunately, 
typical three-dimensional axisymmetric potentials gener- 
ally do not support a global third integral of motion I3, 
even though most orbits are regular. There are various 
ways of different sophistication which can be used to pro- 
ceed. First, the most basic and in principle only wholly 



correct procedure is Schwarzschild's (1979) method of 
modeling a stellar system as superposition of individ- 
ual orbits, which are obtained by numerical integration. 
Second, one may ignore the details of the phase-space 
structure and obtain a global third integral by perturba- 
tion techniques or torus fitting. In this case, ( p^ can be 
generalized to 



/ 



fi{Ri 



aURi)<yURi) 



exp 



LOzJz 
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with the v ertical action J^ ( Binney 1987 , Dehnen fc Bin- 
ney 1996| ). Third, one may use the vertical energy 



E, 



2 z 



$(i?,z)-$(i?,0) 



for /a. This approach is presumably still very good for 
many applications, since Ez is conserved to reasonable 
accuracy for many disk stars (but not for halo stars). 
Finally, one may restrict the analysis to a potential of 
Stackel form, such that a globally valid I3 explicitly ex- 
ists for all orbits (see Bienayme, 1999, for a recent appli- 
cation of this technique). 

6.3. Computer Programs Available 

The computer programs that I have designed in the 
course of this study are written in C-l — h and contain the 
following features. 

1. A general concept for one-dimensional potentials and 
several implementations (power-law, 7-models, Iso- 
chrone) are given. Among other things, this contains 
routines for numerical orbit integration and evaluation 
of Jh, iOR, and w^. 

2. A general concept of warm-disk DFs of the form (|l^ ) 
is given, which allows for integration of general ve- 
locity moments T,v^v2 = / d\» v]lv^f. It also enables 
computation of f{E, L), Ri{E, L), and the derivatives 
of / w.r.t. the parameter functions S' and cr'j^ (useful 
for fitting / to some data). 

3. Code for sampling phase-space points from a DF ei- 
ther pseudo- or quasi-randomly (the latter to reduce 
noise) according to the algorithm given in Section pi 

4. The iterative algorithm of Section |3.2.1 to improve the 
match between moments and targets is implemented. 

These programs are electronically available from me 
upon request. There is also an interface which allows 
one to use a reduced version of these routines from C 
programs. 
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